On the spread of epidemics in a closed heterogeneous population 

Artem S Novozhilov* 
National Institutes of Health, 8600 Rockville Pike, Bethesda, MD 20894, USA 



Abstract 

Heterogeneity is an important property of any population experiencing a disease. Here 
we apply general methods of the theory of heterogeneous populations to the simplest mathe- 
matical models in epidemiology. In particular, an SIR (susccptiblc-infcctive-removed) model 
is formulated and analyzed for different sources of heterogeneity. It is shown that a het- 
erogeneous model can be reduced to a homogeneous model with a nonlinear transmission 
function, which is given in explicit form. The widely used power transmission function is 
deduced from a heterogeneous model with the initial gamma-distribution of the disease pa- 
rameters. Therefore, a mechanistic derivation of the phenomenological model, which mimics 
reality very well, is provided. The equation for the final size of an epidemic for an arbitrary 
initial distribution is found. The implications of population heterogeneity are discussed, 
in particular, it is pointed out that usual moment-closure methods can lead to erroneous 
conclusions if applied for the study of the long-term behavior of the model. 
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Introduction 



Mathematical modeling of epidemics arguably began with the pioneer work of Ross 35| and 
developed considerably ever since (e.g., 0, S])- Especially significant contribution was made by 
Kermack and McKendrick [2^, two students of Ross, who considered the situation of micropara- 
site infection, where contacts between individuals are made according to the law of mass-action, 
all individuals are identical, the population is closed, and the population size is large enough to 
apply a deterministic description (for a brief review see [3] ) ■ Additionally, if it is assumed that 
the individuals are infected for an exponentially distributed period of time, then a usual SIR 
model in the form of ordinary differential equations (ODEs) can be written down for the sizes 
of susceptible, infective and removed classes. 

Since the work of Kermack and McKendrick a great deal of mathematical models were sug- 
gested that relax one or more assumptions that led to the Kermack-McKendrick model. A 
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substantial part of this work was devoted to incorporate heterogeneity into the mathematical 
models, which is also the main subject of the present text. In what follows we retain the as- 
sumptions of random mixing, no inflow of susceptible or infected hosts, exponentially distributed 
infectious period, and validity of deterministic description. Specifically, we will look into het- 
erogeneity in disease parameters (such as susceptibility to a disease); disease parameters are 
considered as an inherent and invariant property of individuals, whereas the parameter val- 
ues can vary between individuals. We analyze heterogeneity that was termed "parametric" by 
Dushoff [8] not addressing important topics of heterogeneity mediated by a structured variable, 
such as explicit space or age structure. 

The most common way to take into account parametric heterogeneity is to divide population 
into groups 0, 16, 17, 3]. An important disadvantage of the subgroup approach is that hetero- 
geneity within a group cannot be incorporated. Another approach, which we al so p ursue, is to 
consider the population as having a continuous distribution (see, e.g., [a, 

or a very large number of subgroups as it was done by May et al. [3l|] (eventually they used a 
continuous gamma-distribution) . 

Our approach to formulate mathematical models is close to that applied in, e.g., where 
known experimental data forced the authors to take into account heterogeneity among hosts in 
their susceptibility to the virus among other key details, and a simple SIR model was adjusted 
to account for new information. The major novelty of the present text is to introduce the well 
developed theory of heterogeneous populations into the epidemiological modeling. Using simple 
models we are able to obtain known results with less effort, and, more importantly, produce new 
analytical results. In particular, we show that any heterogeneous SIR model can be reduced to 
a homogeneous one with a nonlinear transmission function, and present the exact form of this 
function. It turns out that widely applied nonlinear transmission function in the form of power 
relationship, f3S^I'^, can be a consequence of the intrinsic heterogeneity in susceptibility and 
infectivity parameters. For a heterogeneous SIR model the equation of the final epidemic size is 
found. The explicit form of the final epidemic size emphasizes and illustrates the fact that the 
goal to model the evolution of a heterogeneous population for a long time can be accomplished 
only if the exact initial distribution is available. Any moment-closure methods may lead to 
erroneous estimates. 

Our paper is organized as follows. In Section 2 we formulate the basic models and discuss 
various assumptions that might lead to them. In Section 3 we review the necessary analytical 
tools from the theory of heterogeneous populations. In Section 4 a homogeneous model that is 
equivalent to the heterogeneous one is explicitly constructed, and it is shown that the former 
has a nonlinear transmission function; moreover, the well known power transmission function 
is shown to be a consequence of the initial gamma-distribution. In Section 5 the influence of 
heterogeneity on the disease course is studied for an SIR model with distributed susceptibility, 
in particular, the final size of an epidemic is found for an arbitrary initial distribution. Section 
6 devoted to discussion and conclusions. Finally, in Appendix we collect the definitions of the 
probability distributions used throughout the text together with some auxiliary facts from the 
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general theory of heterogeneous models. 



2 The basic models with population heterogeneity 
2.1 The model with distributed susceptibihty 

Suppose that each individual of a (sub-)population has its own value of a certain trait (which 
can be, e.g., susceptibility to a particular disease, social behavior, infectivity, or a hereditary 
attribute) that describes his or her invariant property and has a marked influence on the disease 
course; that is, the key parameters that determine disease evolution depend on the trait values 
and we can speak of the trait distribution or the parameter distributions (in general, we speak 
of a distribution when no ambiguity is expected). The trait value remains unchanged for any 
given individual during the time period we are interested in, but varies from one individual to 
another. Any changes of the mean, variance and other characteristics of the trait distribution 
in the population (or the parameter distributions) are caused only by variation of population 
structure. Such situation is obviously closer to the reality then the usual assumption of the 
population uniformity in the SIR-like models described by ODEs. 

For the moment we assume that the susceptible subpopulation is heterogeneous, and denote 
s{t,uj) the density of susceptible individuals at time t having trait value uj (i.e., the size of 
subpopulation of susceptible hosts having trait values in the range from u; to u; + da; is approx- 
imately equal to s{t,uj)duj, and the total size of the susceptibles S{t) = J^s{t,uj)diJ, where is 
the set of trait values). Assuming that the subpopulation of infectives is homogeneous (under 
the modeling situation), the total size of the population is constant, the disease course keeps 
within the simplest situation "susceptible" — > "infective" — > "removed" , the contact process is de- 
scribed by the so-called mass-action kinetics (i.e., the contact rate is proportional to the total 
population size 0, Q, H), and that the rate of change in the susceptibles is determined by 



the transmission parameter, which is a function of trait values, we can write down the following 
equation for the change in the susceptible subpopulation having trait value lv: 

^s(i,w) = -/?(u;)s(t,a;)/(t), (1) 

where I{t) is the size of the subpopulation of infectives, and P{uj) incorporates information on 
the contact rate and the probability of a successful contact. Hereinafter we assume that the trait 
that characterizes susceptible individuals is the susceptibility to the disease, although it can be 
assumed that different individuals have different contact rates (in the latter case it becomes 
difficult to interpret the equation for the infectives, see below). 

The change in the infective class, if the length of being infective is distributed exponentially 
with the mean time I/7, is given by 

^/(t) = I{t) j (3{u;)s{t, u;) doj - -fl{t) = ^(t)5(t)/(t) - jl{t), (2) 
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where we used the useful notations 

s{t, LO) 



P{t) = / P{uj)ps{t,uj)duj, Ps{t,Uj) 



S{t) 



Hence, /3(t) is the mean value of the function (3{uj), and to has the probability density function 
(pdf) Ps{t,uj) for any time moment t. We need to supplement the model ([I])-© with the initial 
conditions 

s{0,oj) = so{oj) = SopsiO,uj), /(O) = lo, R{0) = 0, (3) 

and with the third equation for the removed dR{t)/dt = ^I{t). Here So and /q are the initial 
sizes of the susceptibles and infectives respectively, and ps{^-,io) is the given initial distribution. 

The model ([I])-® is the basic model we study in this paper. This model was formulated from 
the first principles, as it was done for conceptually similar models in for the transmission 
of virus in gypsy moths, in [33[| for the effect of antimicrobial agents on microbial populations, 
in 



311] for the spread of HIV in the human population, and in [39| for a class of SIS models 



(we note, however, that in the last two examples the frequency-dependent transmission was 
employed [13] , and the heterogeneity of the contact rates was modelled) . 

The model ([II-(l3]) can be also deduced from the general epidemic equation (see (^) 

—s{t,uj) = s{t,uj) / A{T,uj,r]) — s{t-T,r])dTdrj, (4) 

where A(t, u, rf) is the expected infectivity of an individual that was infected r units of time 
ago while having trait value r/ towards a susceptible with trait value uo. If we assume that 
A{t, uj,r]) = I3{uj) exp {— 7t}, and set 

= - / exp{-7T}— s(t - r, rj) drdrj, 

after some algebra we obtain ([l|)-([3]) (see also 0])- As a side remark we note that letting 
A(t, jjj, n) = Piuj)x{T — r), where x{t) is the Heaviside function, we obtain the model studied 
in 0. 



2.2 Model with distributed infectivity 

It can be also assumed that the population of infectives is heterogeneous. Now let /9(u;) be 
the infectivity of an individual with trait value w, and i{t,uj) be the density of the infective 
hosts with trait value uj at time moment t, I{t) = j^i{t,uj) dw. For simplicity we assume 
that the susceptible hosts are homogeneous. The change in the infective subpopulation should 
incorporate the law that specifies which trait value is assigned to a newly infected individual, 
and can be described by the following equation: 

Q fill, 

—i(t,uj)=S(t) / 7p(io,io')p(uj')i(t,uj') duj' — 'yi(t,uj), (5) 
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where ip{uj,uj') is the probabihty that a newly infected individual gets trait value uj if infected 
by an individual with trait value lo' . The change in the population S{t) is given by 

pit) = -msmt), (6) 

where now = j^(3{uj)pi{t^uj) duj^ and pi{t^uj) = i{t,uj)/ I{t). 

The need to specify function ^{lo,u)') precludes the interpretation of the function in 
Section 12.11 as the heterogeneity in the contact rates. It was assumed that the infectives are all 
identical, and thus the supposition that an individual that have had the trait value oj turns into 
another identical infective host is not warranted, at least within the framework of the simple 
model ©-dl. 

A variety of choices for the function ipicOjUj') is possible, but we especially interested in the 
particular case when a newly infected individual gets the same trait value that was possessed 
by the individual who passed the infection, namely 

^p{uJ, to') = 6{uj' — uj), 

where S{uj) is the Dirac delta function (this is similar to the assumptions made in ^9]). In this 
case the equation ([5]) simplifies to the equation, which is very similar to ([1]): 

uj) = f3{Lo)S{t)i{t, u) - ^i{t, to). (7) 

Model Q-dZj) is another example of a simple mathematical model for the spread of an 
infectious disease in a closed population with heterogeneities. The list of possible models can 
be easily extended. For instance, it is straightforward to assume that the parameter 7 is not 
constant for the infected individuals, but rather is distributed with a known initial distribution. 
In this case the equation for the infectives takes the form 

d 

— z(t,a;) = f3S{t)i{t,uj) - j{uj)i{t,u;), 

where /? is now constant. Another obvious generalization is to assume that several model pa- 
rameters are distributed. 

Models ([I])-(l3|) and ©-([Tl) are infinite-dimensional dynamical systems where the evolutionary 
operator specifies complex transformations of the initial distributions. Such models are less 
amenable to qualitative, quantitative, or numerical analysis than their finite-dimensional analogs 
formulated in terms of ODEs. A usual practice is to formulate an infinite-dimensional system 
of ODEs for which some approximations methods (e.g., assuming that the initial distribution is 
close to the normal distribution Q), moment-closure ( 0, Q and B) 

, or numerical methods 



( |39l] ) can be applied. We show below that in some particular cases, when the analytical form 
of the heterogeneous models meets certain requirements, the initial model can be reduced to a 
low-dimensional ODE model, which, in turn, can be effectively analyzed. 
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3 The necessary facts from the theory of heterogeneous popu- 
lations 



To keep the exposition self-contained and for the sake of convenience of references we briefly 



survey the necessary results from |21l . I23l | . We present the results in the form suitable for our 



goal noting that more general cases can be analyzed [23l |. For the proofs we refer to 25|], where 
similar models are considered. Some additional facts are given in Appendix. 

Let us assume that there are two interacting populations whose dynamics depend on trait 
values uji and 102 respectively. The densities are given by ni{t,uji) and n2(i,a;2), and the total 
population sizes Ni{t) = J^^ni{t,uj) duji and N2{t) = J^^n2{t,uj) dui2- Obviously, more than 
two populations can be considered, or some populations may be supposed to be homogeneous; 
we choose two not to be drowned in notations. Assume next that the net reproduction rates of 
the populations have the specific form which is presented below: 

d 

—ni{t,uji) = ni{t,uji)[fi {Ni , N2, ^2{t)) + ^fi {uJi )gi (iVi , iV2 , ^2 (*) )] , 

I' (8) 

— n2(t, wi) = n2(t,a;2)[/2(iVi, A^2, M*)) + M^2)g2{Ni, N2,Mt))], 

where fiiuJi) are given functions, (pi{t) = Lpi{uji)pi{t,u}i) duji are the mean values of 

and pi{t,uji) = ni{t,uJi)/Ni{t) are the corresponding pdfs, i = 1,2. We also assume that ipi{uJi), 

considered as random variables, are independent. The system ([8|) plus the initial conditions 

ni{0,uji) = Ni{0)pi{0,uJi), i = 1,2, (9) 

defines, in general, a complex transformation of densities ni{t,LOi). For the approach to study 
such systems based on the analysis of abstract differential equations in Banach spaces we refer to 
[J) S) S] • Another approach to analyze models in the form ^ was suggested in [2l| . The latter 
is more attractive because eventually one has to deal with systems of ODEs of low dimensions 
(examples of model analysis are given in 

0, 0, MS)- 



Let us denote 

Mi{t, A) = / e^'P^^'^^^p^it, uji) dui, i = l,2, 



in, 

the moment generating functions (mgfs) of the functions ipi{uJi), Mj(0, A) are the mgfs of the 
initial distributions, i = 1,2, which are given. 

Let us introduce auxiliary variables qi{t) as the solutions of the differential equations 

dqi{t)/dt = g,{Nl,N2,^^+l{t)), q^{0) = 0, i = 1,2, (10) 

where indexes that exceed 2 are counted modulo 2. 
The following theorem holds 
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Theorem 1. Suppose that t £ [0, T), where T is the maximal value oft such that has a 

unique solution. Then 

(i) The current means of ipi{uji), i = 1,2, are determined by the formulas 



_,, dMi{0,X) 



dX 

and satisfy the equations 



1 



(11) 



j^^iit) = giiNi,N2,^i+iit))af{t), (12) 

where crf{t) are the current variances of ipi{t,uJi), i = 1,2. 

(ii) The current population sizes Ni(t) and N2{t) satisfy the system 

j^Ni{t) = Ni{t)[fi{Nl,N2,^^+l{t)) + ^^{t)gi{Nl,N2,^^+l{t))], f = 1,2, (13) 

where indexes that exceed 2 are counted modulo 2. 

Theorem [T] gives a method of computation of the main statistical characteristics of (^^(wj); 
the analysis of model ©-([I]) is reduced to analysis of ODE system (fTO|) . (fTT]) . p^ . the only thing 
we need to know is the mgfs of the initial distributions. It is worth noting that the evolution of 
distributions can also be analyzed j23l |. 

Concluding this sections we note that, with obvious notation changes, models ([I])-® and 
©-(IT]) fall into the general framework of the master model ([8]). 



4 Homogeneous models with nonlinear transmission functions 
4.1 Reduction to a system of ODEs 

We start with the equation ([1]), other equations in the system can be quite arbitrary, e.g., the 
full system can contain the class of exposed or several infective classes. Here we show that the 
heterogeneous model that contains ([T|) as a modeling ingredient can be reduced to a homogeneous 
model with a nonlinear transmission function whose explicit form is determined by the initial 
distribution. This result is close to the analysis in [33] where it was argued that the model with 
heterogeneities can be encapsulated in a homogeneous model. Due to the fact that the model 
considered in [s^ is substantially more general than the models we consider here, no explicit 
formulas were found. In the case of model ([l]) nonlinear transmission function can be found in 
the exact form for different initial distributions. 
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According to Theorem [T] we can rewrite equation ([T]) in the form 

d 



dt 



S{t) = -P{t)S{t)I{t), 5(0) = So, 



-q{t) = -I{t), q{0)=0, 



(14) 



dM{0, A) 



dX 



X=q 



where dots denote equations that govern the dynamics of other subpopulations, e.g., these can 
be usual equations for infected and removed classes. M(0, A) is the given mgf of ps{0,uj). 



Proposition 1. Model (I14p is equivalent to the following model: 

j^s{t) = -h{s{t))m, 



(15) 



"dM-i(0,0 




-1 


dC 







where dots denote the same as in ([L 

h{S) = So 

and M^^(0,0 is the inverse function to mgf AI{0, X). 

Proof. The first equation in (I14p can be rewritten in the form 



(16) 



From (dH) p{t) can be represented as (3(t) — ^ 



gives 



dlnM(0,A) 

dlnS{t) d 



X=q(t) 

lnM(0,(?(t)), 



, which, together with the previous. 



dt dt 

or, using the initial conditions S{0) = So, q{0) = 0, 

5(i)/5o = M(0,g(t)), (17) 

which is the first integral to system (jl4p . Knowledge of a first integral allows to reduce the 
order of the system by one. Since M(0, A) is an absolutely monotone function in the case of 
nonnegative P{co) ^ 0, then it follows that 



{t) = M-'iO,Sit)/So) 



(18) 
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where M'^iO, M{0, X)) 
Putting ([H]) into (jT 



A for any A. 

gives 

d _ dM(0,X) 



Sol it), 



A=A/-i(0,S(t)/5o) 
or, by the inverse function theorem, (jlSp with (jl6p . 



□ 



The simple properties of the nonhnear incidence function h{S) are h{0) = 0, /i(S'o) = 
Som, h'{S) > 0. 

Let us consider several examples. Definitions for the probability distributions we use can be 
found in Appendix. In all examples it is assumed that (3{uj) = uj, i.e., the transmission coefficient 
takes the values from the domain of uj with the probability corresponding to Ps{t,uj). 

If the initial distribution is a gamma-distribution with parameters k and ly we obtain 



h{S) = ^ 



s_ 



l/k 



(19) 



If /c = 1 (i.e., the initial distribution is exponential with mean then h{S) = S'^/{vS^ 



h{S) = vS 



(20) 



Expression (jl9p was first obtained in and later used as a nonlinear incidence function in 10(] . 

If the initial distribution is an inverse gaussian (Wald) distribution with parameters and 
V, then 

--Inf- 

In a similar vein other initial distribution can be analyzed, but not all of them have an explicit 
expression for mgf. Other examples of possible initial distributions are uniform distribution 
with parameters a and 6, lognormal distribution, beta-distribution, Weibull distribution, Pareto 
distribution and many others (see Fig. [T]for four different probability distributions). 

The theory we briefly presented in Section [3] is equally valid both for continuous and discrete 
distribution. For instance, if the initial distribution is the Poisson distribution with parameter 
9, then we obtain 



h[S) = S 



(21) 



Note that in the case of the Poisson distribution there is non-zero probability that a randomly 
chosen individual has zero transmission coefficient, i.e., the Poisson distribution implies that 
some individuals are immune. 



4.2 Derivation of the power transmission function 

In standard epidemiological models the incidence rate (the number of new cases in a time unit) 
was frequently used as a bilinear function of infective and susceptible populations: oc SI. In 
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Figure 1: The functions h{S) given by ()16p for four different probability distributions with the 
same means and variances. The diagonal shows the same function for the homogeneous model 
h{S) = pS 



addition to this it is usually argued that there are a variety of reasons that the standard bilinear 
form may require modification, including the assumption of heterogeneous mixing [2^. [sH]. We 
refer to a review paper on the subject 84] for a general account of different models for incidence 
rates, while noting that one of the most widely used models has the form 



(22) 



and is of direct interest to our study. The incidence rate in the form (j22p was first used in 
37[, with the restriction that p, q < 1, but generally it is only required that p,q > 0, see 



also [15|, 127], |29|]. It is interesting to note in the context of our exposition that the exponents 
p, q in (122]) were dubbed as "heterogeneity parameters," but the models itself is considered 
phenomenological and lacking mechanistical derivation [32l |. 

A special case of (122p is g = 1, which was considered in 13, [s^; the values for parameter 
p were considered p > 1. Comparison of the model (|22p with the ODE (|15p when the initial 
distribution of the susceptible subpopulation is a gamma-distribution (see (fT9]) ) let us state the 
following corollary. 



Corollary 1. The power relationship (I22p of the incidence rate in an SIR model for the case 
q = l,p > 1 can be obtained as a consequence of the heterogenous model ([I])-© when the initial 
distribution of the susceptible subpopulation is a gamma-distribution. 
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Let us assume that not only the susceptibles are heterogeneous for some trait that influences 
the disease evolution, but also the infectives are heterogeneous, and consider the simplest possible 
SI model. Let s{t,uji) and i(t,u!2) be the densities of the susceptibles and infectives respectively, 
here we assume that the traits of the two classes are independent, i.e., (3{uji,uj2) = /3i(wi)/32('^2)- 
The number of susceptibles with the trait value uji infected by individuals with trait value 102 is 
given by Pi{uJi)s{t,uJi)P2{^2)i{t,uJ2), and the total change in the infective class with trait value 
^2 is /?2('^2)*(ii ^^2) Jq^ (3ii^i)s{t,u!i) duji; an analogous expression applies to the change in the 
susceptible population. Combining the assumptions we obtain the following model: 



— i{t, UJ2) = P2{^2)i{t,U}2) J f3l{uJl)s{t,UJi)duJl = (52{uJ2)i{t,UJ2)(il{t)S{t). 

Model (f23|) is supplemented with initial conditions s(0,a;i) = SqPs{^-,uji), i(0,u;2) = /oPi(0,W2)- 
In (j23p it is assumed that if an individual having trait value loi was infected by an individual 
with trait value W2 he or she becomes an infective with trait value u)2 (see Section \2.2\\ . The 
global dynamics of (j23p is simple and is similar to the simplest homogeneous SI model. 

According to Theorem [T] the system (I23p can be reduced to a four-dimensional system of 
ODEs. Reasoning exactly as in the proof of Proposition [1] we obtain 

Proposition 2. The model (I23p is equivalent to the model 



where hi{x) is given by (fTHl) . 

Combining Proposition [2] with (jl9p we get 

Corollary 2. The power relationship (j22p of the incidence rate in an SI model for the case 
q > l,p > 1 can be obtained as a consequence of the heterogenous model (|23p when the initial 
distributions of the populations of susceptibles and infectives are gamma- distributions. 

Consequently, it turns out that the power relationship ([2^ . at least for the case p, q ^ 1, 
can be explained on the mechanistic basis by the inherent heterogeneities of the population. Its 
exact form is the consequence of the initial gamma-distributions, but we note that any of the 
transmission functions given in the previous subsection can be well approximated by (jl9p (see 
also Fig. d]). 



^^s{t,u.^) 




(23) 



-S{t) = -h^{S)h2{l) 

^m = h^{s)h2{i\ 
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5 The influence of population heterogeneity on the disease course 



Here we mainly restrict our attention to the model ([I])-® and study its global behavior. First 
we state almost obvious proposition: 

Proposition 3. Let Si{t) = J^si{t,uj) duj be the solution of ([I])-® with the initial condition 
si{0,uj) = SqPi{0,uj), and 5*2 (t) = f^S2{t,uj) dto be the solution of (H])-® with the initial 
condition S2{0,lu) = SoP2{0,uj), such that /3i(0) = /3i(0) and (Tf(0) > cr^(O), where /3i(0) = 
f3{uj)pi{0, u) dio and af{0) = f^{(3{uj) — /3j(0))^pj(0, w) dw, i = 1,2. Then there exists e > 
such that Si{t) > 52(t) for all t £ (0, e). 

The gist of this proposition is very simple: the more heterogeneous the susceptible hosts the 
less severe the disease progression under the model ([I])-®. 

Proof. Differentiating the first equation in (I14|) and using (I12p we get 

S"{t) = I^Sia\t)+Pit))-p{t)rS, 

or, at the initial time moment, S'"(0) > 5*2 (0). Since S{t) is continuous the proposition follows. 

□ 



We remark that this proposition also holds for more general model (I14p . Moreover, we can 
replace equation ^ with equation 

^s(t, io) = -(3{io)s{t, oj)I{t) + s{t, uj)[. . .] 

where dots denote terms describing demography, migration or the lost of immunity by removed 
individuals, the only condition is that these terms cannot depend on Even in this case 

Proposition [3] still holds. For the model ([S])-® the opposite proposition is true: the more 
heterogeneous the infective class in infectivity, the more severe the disease progression, which 
follows from the fact that i3'{t) > in this case, and, consequently, 5*1 (0) < 5*2(0). 

It is interesting to note that in the case of proportional mixing (frequency-dependent trans- 
mission) knowledge of only the initial variances of the parameter distributions does not allow 
inference on the short ran behavior [sO^. 

One of the main characteristic of SIR models is the final size of the disease, which is often 
expressed in the number (or proportion) of susceptibles that never get infected. For the Kermack- 
McKendrick model 

dS/dt = -I3SI, dl/dt = PSI - 7I, dR/dt = 7/, 

it is well known that the desired number, which we denote as ^(oo), is given by the root of the 
equation 

5(00) =5oexp|-^(A^-5(oo))| (24) 



12 



on the interval (0,50). Here N is the constant size of the population. It is easy to show that 
this root always exists. 

Recall that M(0, A) is the mgf of the initial parameter distribution. For the model ([I])-® 
the following theorem holds. 

Theorem 2. The size of the susceptible subpopulation that escapes infection within the frame- 
work of the model is given by the solution of the equation 

S{oo) = SqM (0, {S{oo) - N)j-^) (25) 

satisfying condition < S'(oo) < 5*0. 

Proof. First we note that exactly as it was done in Proposition [H we can reduce the system 
([I])-© to the system 

dS{t)/dt = -h{S{t))I{t), 

dl{t)/dt = h{S{t))I{t) - jl{t), (26) 
dR{t)/dt = "/I{t). 

Using (I16p and dividing the first equation in ()26p by the third one we obtain 



dM-i(0,e) 



"7 du. 

-S/So 



So 



dC 

Integrating from to oo gives 

Using the identities R{oo) = N — S{oo) (since /(oo) = 0) and M^^{0, 1) = we obtain 

M-(S(oo)/S,) = -^^^^, 

from which (j25p follows. Due to the fact that M(0, A) is an increasing function, the solution of 
(j25]) satisfying < S{oo) < Sq is unique. □ 



Remark. If we consider nondistributed parameter (formally, we can let P{uj) = (3 = const, or, 
equivalently, s(0,u;) = Sq5{uj — w), /3((D) = const, where 5{ijj) is the delta- function) , we obtain 
(1^ from (125]). 

Arguing in the same spirit as it is done in 0] (e.g., p. 183), the problem of the epidemic 
invasion can be considered. Assuming that initially all the population is susceptible (formally. 
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for our model, S{—oo) = N, q{—oo) = 0), from ([25]) the equation for the fraction of susceptible 
population that does not get infected follows: 



z = M(0,-N{1 - z)/7). 



(27) 



Here z = S{oo)/N. Equation (j27|) always has the root z = 1. If the basic reproductive number 
0], defined here as 



satisfies the condition Rq > 1, then there is another root of (j27p in the interval < z < 1. 
This root gives the sought fraction. The proof of the existence of this root under the threshold 
condition iig > 1 is straightforward and can be conducted similar to the homogeneous case 
(e.g., [5]). This result can be illustrated by the case when the initial distribution is exponential 
with parameter u: the equation for z is quadratic: Rqz^ — {1 + Rq)z+1 = 0, where Rq = N/i^ju). 
This equation has the roots 1 and I/Rq- If i?o > 1 then the fraction of susceptible population 
that escapes the disease is I/Rq. 

Comparing the results obtained for the heterogeneous SIR model ([I])-© with the well known 
results for the simple homogeneous SIR model, we can conclude that the questions of the disease 
invasion can be studied in the framework of the homogeneous model because population hetero- 
geneity does not impact the basic reproductive number ([28]) (this holds, obviously, if we identify 
the mean value of /3(w) over the population of susceptibles at the initial time moment with the 
usual constant /5 in the homogeneous model). From the other hand, the heterogeneity of the 
population has direct impact on the final size of the disease since equation ()27[) depends on the 
initial distribution in contrast to the homogeneous analogue z = exp{— iio(l ~ z)} (surprisingly, 
the last formula is valid under variety of different conditions [30]). 

In Fig. [2] the final size of the susceptible population versus the initial variance of the pa- 
rameter distribution is shown for four different initial distributions that have the same means 
at t = 0. From Fig. [2] it can be seen that the more heterogeneous the population of suscep- 
tibles, the less severe the disease not only in the short run (Proposition [3|) but also globally 



(see [SjJ for the same result for frequency-dependent transmission). At the same time, it is 
worth emphasizing that the conditions /3i(0) = P2{0), crf{0) > cjKO) for two different initial 
distributions do not imply that zi > Z2, where Zj are the solutions of ()27p . A counterexample 
can be easily found (e.g., taking gamma-distribution with parameters A; = 2, = 4 and uniform 
distribution on the interval [0,1], N/'j = 20, we find that zi = 0.093 < Z2 = 0.112 whereas 



If we compare two distributions of the same family then sometimes it is possible to prove 
rigorously that, on the assumption of equal initial means and different second moments, more 
heterogeneous population (i.e., the one that has larger initial variance) experiences less severe 
disease. This is true, e.g., for two gamma-distributions: 




ai{0) = 1/8 > (T2(0) = 1/12). 
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Figure 2: The size of the susceptible population that never gets infected, 5(oo) versus the initial 
variance of the parameter distribution ct^(O) for four different initial distributions, ,5(0) is the 
same in all cases. The parameters are S{0) = 999, 1(0) = 1, 7 = 20, /3(0) = 0.05 



Proposition 4. Assume that we have model (HI) -([3]) with two initial gamma-distributions with 
parameters ki,ui and k2,i^2 such that I3i{0) = (32{0) and crf{0) > (t|(0). Assume that Rq > 1. 
Then solutions of (|27|) that belong to (0, 1) are always satisfy z\> Z2. 



Proof. We need to prove that 



{k2 + Roil - z))"^ 



(fcl + Ro{l - z))'=i k''' 



> 1 



for any z € (0, 1). This follows from the fact that /(r) 
increasing function for r > 0, i.e., f'{r) > 0. 



(^2 + r)^^ /{ki + r)^^ is monotonically 



□ 



Remark. The last proposition can be extended to other initial distributions, e.g., it holds for 
Wald distribution, for a uniform distribution and some others. However, for any distribution 
that is not determined by its first two moments (i.e., that depends on more than two parameters) 
an analogous proposition is no longer true. 



6 Discussion and Conclusions 

Here we have presented several results concerning the course of a disease in a closed heteroge- 
neous population, where heterogeneity is mediated by invariable traits whose distributions are 
determined by population structure at every time moment. One of the purposes of the present 
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text is to introduce in the area of epidemiological modeling the general technics of the theory 



of heterogeneous populations [2l|, |23| which allows reduction of the initial infinite-dimensional 
dynamical system to an ODE system of a low dimension. 

A usual strategy in the literature when analyzing systems similar to ([I|)-(l3]) is to consider an 
infinite-dimensional system of ODEs where the variables can be moments or cumulants of the 



corresponding distributions (e.g., 111. l33l. l39l|) and then analyze this system, or consider one of 
many possible moment-closure strategies to extract valuable information. Having the general 
theory outlined in Section [3] and the known results from the literature, we can critically discuss 
the latter and be prepared to possible pitfalls. 

For example, in [ll]] the equation for the final epidemic size was obtained by using two 
different strategies: first, an initial gamma-distribution was assumed and analytical treatment 
was applied, and second, using the procedure suggested in an approximation method was 
used in which the infinite-dimensional system of ODEs was replaced with two equations under 
the assumption that the coefficient of variation is constant. It is not surprising that Dwyer 
et al. obtained identical results because the gamma-distribution, according to the theory of 
heterogeneous populations, is the only continuous distribution which does not change its shape 
during the system evolution, and keeps the coefficient of variation constant (see Appendix). 
Therefore, the conclusion that "...the assumption of gamma-distributed susceptibility is not 
strictly necessary to derive equation [for the final epidemic size]" is not valid in many situations. 
Another initial distribution, how it is shown by ()25l) . can yield another equation for the final 
epidemic size and, consequently, can produce significant discrepancy with the moment-closure 
approximation suggested by Dushoff [sj. 

The well known fact from the theory of heterogeneous populations that to model the system 
dynamics for a substantial time period we need to know the exact initial distribution implies that 
any results obtained for an epidemic in a heterogeneous population on the ground of knowledge 
of only several first moments of the initial distribution have to taken with extreme care. One, 
two, or more first moments of the initial distribution can be insufficient or even misleading. We 
also note that the short run behavior can be predicted when we have information only on several 
first moments (Section [5]). 



Another initial distribution which was used in the literature is the normal distribution |33l |. 
For the normal distribution we have that = 0, i ^ 3, where Kj is the i-th cumulant. Combin- 
ing the last property with the fact that the initial normal distribution remains normal within 
the framework of heterogeneous models we obtain that Ki{t) = 0, i ^ 3, holds for any t (see 
Appendix). This was used in [33] to obtain an explicit solution to the equation 

d 

—n{t, r) = Kgn{t, r) — rn{t, r). 

Here n(t, r) is the number of cells at time moment t, which are killed by antimicrobial agents 
with the kill rate r, and Kg is a constant (notations are changed from the original). First we 
note that, using Theorem (H we obtain explicit solution of this equation for an arbitrary initial 
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distribution of the kill rate: 

N{t) = NoM{0, -t) e^p{Kgt}, 

where N{t) = n{t, r) dr. Second, the results from Sections H] and [5] cannot be applied to 
the normal distribution because this distribution is defined from — oo to oo and thus the cor- 
responding mgf does not have an inverse. Which is more important, however, the total system 
dynamics can be influenced by these negative kill rates even if they occur with vanishingly small 
probability (we note that this issue is discussed in [l^). An example of such influence can be 
found in 2j] where infinitely large growth rates occurring with small probabilities drive the 
population to explosion. Therefore any approximations based on an initial normal distribution 
in a situation where parameter can take only nonnegative values should be taken with care. 

Summarizing the main results we can assert that the theory of heterogeneous populations 
can be successfully applied to many different mathematical models in epidemiology. Examples 
are given in Section [3j In many simple cases the original model can be reduced to a model 
described by ODEs, which simplifies the analysis. The law of mass action for a distributed 
susceptibility model implies a nonlinear incidence function in a homogeneous model. Moreover, 
one of the well known transmission functions, power relationship, follows in exact form from 
the initial gamma-distribution, at least in the case when exponents exceed one (see Section S]). 
Therefore, a mechanistic derivation has been given to the transmission power function, which 
was shown previously approximate real data with high accuracy. The short term behavior of the 
models considered can be approximately described knowing only two first moments of the initial 
distribution, whereas the long-term behavior depends on the exact initial distribution and can 
vary significantly (Section [5]) even for the distributions whose several first moments are identical 

(Fig. m- 

It is a tempting challenge to include various demography processes to the analyzed models. 
The main obstacle is the need to specify the function </>(ti;, uj') similar to the one used in The 
delta-function yields the models that can be analyzed using the general approach from Section 
[3] but it is usually difficult to interpret the underlying assumptions. These problems are the 
subject of ongoing research. 



A Appendix 

In Appendix we collect the definitions of the distributions used throughout the main text. The 
definitions are taken from [l^ and In addition to that we list some facts concerning 

evolution of these distributions if they are used as the initial distributions for the models studies 
in the text. Everywhere below it is assumed that /9(w) = uj. Inasmuch as we are interested in 
characteristics of distributions depending on time, the following formula is very useful (see [2li|): 



M(M+^ (29) 
^' ^ M{0,q{t)) ^ ' 
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where q{t) in the solution of the corresponding auxihary differential equation (see Theorem [T]) , 
M(t, A) is the mgf of the parameter distribution at time t. Equation ()29p shows that the mgf at 
any time instant can be expressed using the initial mgf. 
Gamma- distribution 

The pdf of gamma-distribution with parameters k and z/ is given by 



m 



UJ 



k—l„—uu} 



UJ ^ 0, k > 0, 1/ > 0. 



(30) 



The mgf of gamma-distribution is M{0, A) = (1 — X/v)~^ . 

It follows from ()29p that for i > the distribution does not change its form, i.e., it is gamma- 
distribution with parameters k and s — q{t). The mean and variance of the distribution are given 
by 

k 2/ 



Pit) 



a\t) 



k 



s-q{t)' {s-q{t)y 

Note that at any time moment the coefficient of variation is constant: cv = a{t)/f3{t) = Xj^fk. 
Actually, gamma-distribution is the only continuous distribution whose coefficient of variation 
remains constant with time within the framework of heterogeneous models. 
Wold distribution 

The pdf of inverse gaussian (Wald) distribution with parameters and is 



PiO,uj) 



27ra;3 



1/2 



exp 



2fj.^uj 



(uJ - fl)' 



The mgf of Wald distribution is given by 



M(0, A) = exp <^ - (l 



U! > 0, fj,, v > 0. 



1/2N 



Again the distribution remains Wald distribution with parameters 

n 1 

L 

and temporal characteristics of the distribution are 

2fi^q{t) 



Pit) = ^ 
Beta- distribution 



a\t) 



PHt) 



cv 



Pit) 



(31) 
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Sometimes it is useful to study evolution of distribution with compact support. A good 
candidate in this case is the family of beta-distributions with pdf 



1 



a ^ uj ^b, n > 0, r2 > 



(32) 



(6 - ay^+''^- 



where B{ri,r2) is the beta-function. 
The initial mean and variance are 



a\0) 



rir2 



ri + r2 



(n + r2)2(ri + r2 + 1) 



Unfortunately in the case of beta-distribution it is impossible to write down the mgf, and, 
correspondingly, the temporal characteristics of the distribution. In a special case ri = r2 = 1 
we have a uniform distribution on [a, b]. Equation (j29p shows that in this case for t > the 
distribution is no longer uniform but turns into truncated exponential distribution. 

In the text we used beta-distribution on [0, 1]. 

Log-normal distribution 

The pdf of log-normal distribution defined on the nonnegative half-axis uj ^ with parame- 
ters n and u is 



As in the case of beta-distribution we cannot present explicit formulas for the mgf and other 
characteristics. We note that this distribution can be used only if q{t) < for t > 0, otherwise 
the integral in the mgf diverges. This is the case, e.g., for the model (HJ-dS]), but not the case 
for ©-([T]), for which the log-normal distribution cannot be used. 
Normal distribution 




(33) 



The initial characteristics are 



/3(0) = exp{^ + ^^2} 



a^{0) = exp{2// + u'^}{exp{u^} - 1). 



The pdf is 




(34) 



The mgf is M(0, A) = exp {iA(2/i + Ao"^)}. The temporal characteristics are 

Pit) = fi + 2q{t)a^ , a\t) = a^. 



Poisson distribution 
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In the same spirit discrete distributions can be managed. Consider, e.g., the Poisson distri- 
bution with parameter 9, i.e, 

p(0,a;) = — e", 6* > 0, a; = 0, 1, . . . (35) 

then M(0, A) = exp {e{e^ - 1)}, and 

(3{t) = a{t) = i9exp{g(i)}, cv = 1. 

Other possible initial distribution can be considered in a similar vein. 

Acknowledgments. The author thanks Dr. A. Bratus' and Dr. G. Karev for insightful 
discussions and helpful suggestions. 
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